Amplitude coda of classical waves in disordered media 
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The propagation of classical waves in the presence of a disordered medium is studied. We consider 
wave pulses containing a broad range of frequencies in terms of the configurationally averaged Green 
function of the system. Damped oscillations in the time-dependent response trailing behind the 
direct arrival of the pulse (coda) are predicted, the periods of which are governed by the density of 
scatterers. 
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Wave propagation in disordered media is a broad research topic, with many applications. Electron transport 
in mesoscopic systems, light diffusion in opaque media, or acoustic propagation in the subsurface of the earth are 
examples of the interdisciplinary character of this general subject. All these examples share a number of common 
properties, the most important being that they are governed by wave equations. Methods and techniques used in one 
field of study can in principle be applied to others, and such an approach has proved successful in the past. Anderson 
localization (!]] for instance, was originally discovered in the quantum realm of electrons but later understood to be a 
general wave phenomenon 

Other contributions from traditionally quantum methods into the study of classical waves are the diagrammatic 
perturbation technique || and the random matrix theory Q . Whereas the latter explains some universal features in 
the response of a disordered environment which are independent of the detailed structure of the system, the microscopic 
nature of the former allows a more specific modelling of the medium. If one wishes to image the medium by observing 
the signal which is generated by a source and scattered by the system inhomogeneities, the desired information must 
be sought in the non-universal part of the measurements. 

In addition to studies of classical wave propagation based on nearly monochromatic sources, the time-dependent 
response of a short pulse containing a broad range of frequencies is also of considerable interest. This is typically the 
case in seismic studies of the earth. Excited either by earthquakes or by artificially controlled explosions, seismic waves 
travel long distances before being registered by an array of detectors. The measurements consist of time-dependent 
functions describing how the subsurface responds to the excitations. The portion of the wave lagging behind the first 
arrivals (the so-called coda) is delayed by being repeatedly scattered by interfaces and inhomogeneities and therefore 
carries information about the medium |Q . Only recently it has been realized that the coda of seismic waves can be 
used for mapping the subsurface of the earth as well as locating earthquake centres ||] . 

Unlike their quantum counterparts, the amplitude of classical waves can be probed quite trivially. This is the case 
for seismic and acoustic waves in general. Ultrasonic wave amplitudes are also observable and have been used to 
investigate multiple scattering effects in disordered media |3|,[l0| • Amplitude and phase of electromagnetic waves have 
also been measured in the microwave region and used to treat dynamical aspects of the propagation |jll],[l2) . Such a 
wide applicability motivates us to study the time evolution of the wave amplitude within a disordered medium. 

Bearing in mind the microscopic approach mentioned above, we have recently investigated the time response of a 
disordered medium composed of small spheres (Rayleigh scatterers) . The scatterer is characterized by a set of 
resonances which appear as poles of the scattering matrix on the lower half of the complex frequency-plane. A simple 
expression for the time-dependent wave amplitude is obtained when all but the lowest order poles are discarded. We 
subsequently considered an ensemble of randomly placed identical scatterers and related some features of the response 
to the microscopic details of the structure. 

In this paper we extend this idea by accounting for all poles of the scattering matrix. On the one hand we loose 
the simplicity of the single pole approximation, but on the other the single scatterer contribution is better described 
by considering scattering events previously neglected. In spite of the somewhat more complicated formalism, we are 
still able to relate the delayed-time response to the microscopic parameters of the system in the limit of impenetrable 
scatterers. Oscillations in the time-dependent amplitude of the wave are identified, the periods of which depend on 
the concentration of scatterers. 

We consider a general scalar wave equation in the frequency-domain given by the following eigenvalue equation 

{-V 2 + V(r,Lj)}V m (r,u) =E m (u)V(r,u) , (1) 

where \E , rn (f, u>) is the wave field for a frequency uj at position r. Both the energy E m and the potential energy V 

are frequency-dependent and given by V(r, uS) = ^ 1 1 57^ ) and E m (to) — The wave velocity c(r) is position- 

dependent, cq being the velocity for the homogeneous background medium. Associated with Eq. (yj) there is a Green 
function G given by 

G(r,r ,u) = }^— 2 — , (2) 

to % _ E m(v) + ir\ sgn[u) 

where r\ is a positive infinitesimal number. G{r, fo, u>) represents the response measured at position r due to a 
stationary perturbation of frequency u) produced at r*o ■ For the sake of simplicity we assume r*o at the origin throughout 
this paper. The Fourier transform to time-domain gives the physically relevant quantity G(r, t), namely the response to 
a pulse containing the full spectrum of frequencies. From this quantity a simple convolution gives the wave amplitude 
for a pulse of arbitrary shape. 

For a disordered system with many scatterers we focus on the configurationally averaged Green function (G(r, t)). 
In wave vector k and frequency space, it is given by 
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(G(k, u)) = -—2 , (3) 



where the so-called self-energy E(fc,w) is a complex quantity. The solution to the disordered problem is determined 
by the behaviour of the self-energy. For spherically symmetric scattering cross sections (s-wave) and within the 
independent scatterer approximation fl4|| , valid for a low concentration of small defects n, the fc-dependence of the 
self-energy disappears and £ becomes 

= »co(5o(a,)-l) 

where So(w) is the s-wave scattering matrix for the individual scatterers. Note that equations (||) and (|J) bridge the 
gap between the microscopic scale, in terms of So of a single impurity, and the macroscopic medium, described by 
the average Green function (G) . 

The scatterers are modelled by small spherical regions within which the waves propagate with velocity c. The 
corresponding S-matrix is given by fig] 

SoM - -e-*" d/c ° ( 1 + 7~~~~7"\ | • I ) 

V fcot(^)- 

where d is the diameter of the scatterer. So has an infinite number of poles in the frequency-domain, each one of 



them on the lower half of the complex frequency-plane. As shown in Ref. |13|, a simple expression for the response 
is obtained when all but the lowest order poles are neglected. Here we want to include all poles of the S'-matrix and 
investigate the effect this may bring to the final response. After substituting Eq. (|^) into (Q), the Fourier transform 
of Eq. (g) yields the average Green function (G(f, t)) in position and time domain. Alternatively, we can look at the 
scattered response (AG(r,t)), defined as the difference in amplitude between the total signal and the impurity-free 
response, which describes how the scattered portion of the waves evolves in time. Figure |l| shows the calculated 
(AG(r, ()) for a random medium with identical scatterers at n = 5 x I0~ 2 and n — 10 _1 volume concentrations and 
c = and c = g. n and c are expressed in units of (l/(i 3 ) and Co, respectively, whereas the time is in units of d/co- 
The distance r to the source is arbitrarily fixed at r — lOd. As expected, the direct arrival of the wave is followed 
by an exponentially decaying response whose features depend on n as well as on the scattering strength. In addition, 
clear oscillations in the scattered response can be seen in the insets of figure [fl Although not shown in the figure, 
oscillations are also present when c > cq. 

To understand the origin of those oscillations it is helpful to simplify the problem and assume that the wave velocity 
vanishes inside the scatterers, i.e., c = 0. This is equivalent to a set of obstacles which are impenetrable to the waves, 
modelling for instance strong discontinuities in the constitutive parameters. One example is a solid porous medium 
filled with rarefied air, or any other environment with scatterers of relative low sound velocity. Alternatively, in 
the case of electromagnetic waves, obstacles with high refractive index are required. With this simplification, the 
concentration n and the diameter d are the only variable parameters and the S-matrix becomes So = — e ~ tujd / c '> . The 
response (G(k,t)) in wave-number and time domain is an intermediate step towards (G(f,t)) and results from the 
time Fourier transform of Eq. (^), given by 

(G(k,t)) = ±l <!--•[-){- , : }■ «h 
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Although the spatial part of the Fourier transform remains to be done, Eq. is useful for understanding the 
time dependence of the response. The frequency integral can be evaluated by contour integration. The poles are 
at fc-dependent frequency values which govern the actual dispersion relation of the wave motion. In the absence of 
scatterers (n — 0) for instance, the poles are at w = cok and a further Fourier transform into space-domain gives 
the usual free-space solution [[16| Q(r,t) = j^S(cot — r). For the poles acquire a concentration dependence 

u>j(k,n), where j is a labeling index. In general, the residue 1Z associated with a first-order pole u>j(k,n) is given by 

^(k,n) (*■")« 

K^(k,n)} = ^-^ , (7) 
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where <^ stands for the w derivative of <f> = - \ k 2 + g ^ e -^d/c _|_ ij evaluated at w = LJj(k,n). The 

time dependence of the residue is entirely in the exponential e - lu j( fe >™)* . The integral of Eq. (||), written as a sum of 
residues, is 

<G(M)> = __L^ ^*."J« . (8) 

In this equation each term gives an exponentially damped oscillatory contribution, whose period and rate of decay 
depend on the real and imaginary parts of u>j{k,n) 1 respectively. Whether this time dependence persists in the final 
response depends on the spatial Fourier transform. 

After integrating over the angular variables, the remaining transform becomes 

/+oo 
dkA(k,r) e - lu) > (k > n)t , (9) 
j 

where A(k, r) = lk ""^l^^'"" 1 ■ Except for the damping component induced by the imaginary part of the pole, 

the integrand in Eq. (^|) is a periodic function of t, oscillating with a frequency given by the real part of ujj(k,n) 
and with an amplitude governed by A(k, r). For large values of t the exponential oscillates rapidly as a function of 
k and the dominant contribution to the integral comes from the vicinity of /c-points at which uij(k,n) is stationary. 
In other words, the values of k that effectively contribute to the integral are those satisfying the condition = 0. 
Although each pole in Eq. (||) contributes with a different term in (g), those meeting the stationary phase condition 
will dominate the signal. 

The existence of such stationary points is implied by the oscillations found in Fig. |l| and we now try to locate 
them. A closer look at Eq. (Jj|) shows that the poles u>j(k, n) are functions of k 2 . Therefore, k = always satisfies the 

stationary phase condition —r& = 0. We assume this value of k as a possible candidate for explaining the observed 
results in Fig. |l]. The corresponding poles are the solutions of a transcendental equation, which can be further 
simplified if we expand the S- matrix for small u>. Assuming that n is small, which is the basic condition for the 
validity of the independent scatterer approximation, the zero order term is sufficient and we find for the real part of 
the pole that $t[u}j(k — 0,n)] — (n) 1 / 3 . Associated with this pole is a simple expression for the oscillation period 
T given by 

T = 2 7 r/5R[u; i (fc = 0,n)} = ^± . (10) 

The predicted oscillation periods based on the stationary phase argument, shown as a continuous line in figure 
agree very well with the observed periods obtained by numerical integration of the Green function, represented by 
the points. Although the stationary phase approximation provides no information about the respective weight of the 
oscillations, it is sufficient to explain their existence in the time response (G(r,t)}. With such a simple expression for 
the period, the oscillations are easily related to a microscopic structural parameter of the system, namely the average 
separation a = \j \fn between scatterers. 

A simple physical picture for the oscillations in the time-dependent amplitude is as follows. When scattered 
only once, the S'-matrix Sq induces a phase shift in the wave corresponding to a sign change. A second scattering 
event restores the original sign, which is again changed by a third scatterer, and so on. In this way, positive and 
negative contributions to the scattered response occur, depending on the number of scattering events. The maximum 
probability of being scattered a certain number of times depends on the time lapse, and therefore dictates the overall 
sign of the final response. For that reason the sign of the response changes as a function of time and we can also 
understand why the average distance between scatterers appears as the period of the oscillations. On the opposite 
side of the velocity spectrum, infinitely large values of c lead to So = e - lujd / c ° ^ causing no sign change in individual 
scattering events. In this case, no oscillation should occur and, indeed, are not observed in the numerical solution for 
this limiting case. 

Therefore, according to Eq. (JTo|) , it might be possible to obtain the values of a and n, once such oscillations in 
(G(r,t)) are observed. Whether this structural information can be extracted from an actual measurement depends 
on the type of disorder in the system. Within the conditions for the validity of the present model, i.e., tod <C 1 and 
a/d^> 1, the oscillations do not depend on the distribution of the size of the scatterers. 

Well defined oscillations in the time-delayed response have been recently observed in ultrasonic waves transmitted 
across a 2-dimensional disordered structure [Eo[. They have a different origin, but we suspect that the superimposed 
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beating pattern reflects the physics discussed here. This hypothesis could be tested by systematic experiments with 
different concentration of scatterers. 

In summary, we have observed oscillations in the time dependent response of a system composed of randomly placed 
spherical scatterers. Based on the stationary phase approximation, we have identified the origin of the oscillations 
and derived a simple expression for the periods in the limit of impenetrable scatterers, whose dependence on the mean 
separation between obstacles was established. It is suggested that such oscillations should be observable in disordered 
systems providing an indirect way of measuring the density of scatterers. When the scatterers are penetrable the 
simple analytical result presented here has to be modified, but can still be useful in the comparison to the numerical 
results. 
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FIG. 1. Scattered response (AG(r, t)) (in units of 10 x d /co) for r = 10 d. The graphs at the top correspond to the case of 
impenetrable spheres where c — 0. The bottom graphs are for c = co/5. The graphs on the left and right are for n — 5 x 10~ 2 (in 
units of 1/d 3 ) and n — 10 _1 (in units of 1/d 3 ), respectively. The insets are amplified by a factor 10. 



FIG. 2. Oscillation periods (in units of d/co) as a function of the concentration n (in units of 1/d 3 ). The continuous line 
results from the stationary phase aproximation whereas the points are the periods obtained after numerical integration. 
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